Diagnostic imaging is an important tool for observing everything from subcellular structures to the environment, and tomographic reconstruction is an imaging technique that underlies nearly all key diagnostic imaging modalities. Such tomographic modes include X-ray Computed Tomography (CT), Positron Emission Tomography (PET), and Single Photon Emission Count Tomography (SPECT). In addition, certain acquisition methods for Magnetic Resonance Imaging (MRI) use tomography, and new techniques are emerging such as electrical impedance tomography (EIT) and optical tomography. Tomographic reconstruction is also fundamental to numerous other applications in science and engineering, including electron microscopy to determine subcellular structure, nondestructive evaluation (NDE) in manufacturing, geophysical exploration, environmental monitoring, and remote sensing.
An important component in the various tomographic reconstruction procedures is the operation of computing a set of projections of a given image, called forward projection or reprojection. Reprojection is a process by which projections are produced from an image, such that, if the projections are filtered and back projected, they yield the original image. While reprojection is widely used for various purposes, known algorithms that employ the reprojection processes incur a high computational expense.
Reprojection is of interest in several applications. Reprojection is used in X-ray CT, see, C. R. Crawford, J. G. Colsher, N. J. Pelc, and A. H. R. Lonn, "High Speed Reprojection and its Application," Proc. SPIE--Int. Soc. Opt. Eng. Conf. Medical Imaging II, Newport Beach, Calif., vol 914, pt A, pp. 311-18, 1988; in iterative beam-hardening correction algorithms, see, P. M. Joseph and R. D. Spital, "A Method for Correcting Bone Induced Artifacts in Computed Tomography Scanners," JCAT vol 2 No. 1, pp. 100-108, 1978, and U.S. Pat. Nos. 4,217,641 to Naparstek and 5,438,602 to Crawford, et al.); in streak suppression algorithms, see, G. Henrich, "A Simple Computational Method for Reducing Streak Artifacts in CT Images," Computerized Tomography, vol. 4, pp. 67-71, 1980 and U.S. Pat. No. 5,229,934 to Mattson et al.; in algorithms for the removal of artifacts caused by the presence of metallic implants in a patient, see, G. Glover and N. J. Pelc, "An Algorithm for the Reduction of Metal Clip Artifacts in CT Reconstructions." Medical Physics, vol. 8 No. 6, pp. 799-807, 1981 and U.S. Pat. No. 5,243,664 to Tuy; or other high density objects, see, U.S. Pat. No. 4,709,333 to Crawford; and in correcting for missing data, see, U.S. Pat. No. 5,396,528 to Hu et al., and partial volume effects, see, U.S. Pat No. 5,727,041 to Hsieh.
In PET and SPECT imaging, reprojection has been used to compensate for attenuation, see, J. M. Ollinger, "Reconstruction-Reprojection Processing of Transmission Scans and the Variance of PET Images," IEEE Trans. Nuc. Sci. Vol. 39 No. 4, p 1122 August 1992. Reprojection is also used in detection and compensation for various acquisition errors, see, S. C. Huang and D. C. Yu, "Capability Evaluation of a Sinogram Error Detection and Correction Method in Computed Tomography," IEEE Trans. Nuc. Sci. Vol. 39 No. 4, pp. 1106-1110 August 1992, and E. C. Frey, Z. W. Ju, and B. M. W. Tsui, "A Fast Projector-Backprojector Pair Modeling the Asymmetric, Spatially Varying Scatter Response Function for Scatter Compensation in SPECT Imaging," IEEE Trans. Nic. Sci, vol. 40, No. 4, pp. 1192-7, August 1993, including acquisition errors caused by patient motion, see, U.S. Pat. Nos. 4,858,128 to Nowak, 5,552,605 to Arata, 5,579,358 to Lin, and 5,848,114 to Kawai et al.; and in accounting for Poisson noise statistics, see, J. M. Ollinger, "Iterative Reconstruction-Reprojection and the Expectation-Maximization Algorithm," IEEE Trans. Med. Imag. Vol. 9 No. 1, p. 94 March 1990. Reprojection is also a critical component in iterative tomographic reconstruction algorithms, which are the preferred method in imaging modalities such as PET, SPECT, and nondestructive testing. Iterative tomographic reconstruction algorithms are also critical in X-ray CT algorithms for handling a variety of the aforementioned artifacts caused by the presence of metal clips in CT images.
There is a need for fast reconstruction of such tomographic images, especially with the advent of new technologies that are capable of collecting large quantities of data in real time, for example, multi-line spiral CT, cardiac imaging using X-ray CT, and an upcoming CT fluoroscopy. In addition, there exists an increasing demand for real-time interventional imaging, e.g., to monitor and guide surgery. Very often with known techniques the reconstruction of images from such data becomes a bottleneck. While it has been a subject of continuing efforts in industry and academia since the introduction of CT, known efforts to speed up tomographic reconstruction have been unsuccessful. In addition, the development of methods to speed up tomographic reconstruction has been recently identified as a major thrust area by groups such as the National Institute of Health (NIH).
Current iterative methods are far more expensive than the conventional direct filtered back projection (FBP) reconstruction, and typically require orders of magnitude more computation. This has motivated intensive work on the acceleration of iterative methods dating back to the first CT scanner around 1971, which used iterative reconstruction. More recently, the need to develop fast iterative reconstruction methods has been identified as a pressing need, see, Imaging Sciences Working Group, "Matching Clinical and Biological Needs with Emerging Imaging Technologies," tech. rep., Diagnostic Imaging Program, National Cancer Institute, 1998. Because of the key role played by the process of reprojection in the various recontruction algorithms, the importance of accelerating it has been recognized as early as 1978.
Prior methods for reprojection are problematic. Direct methods work directly on the data, and involve direct computation of weighted sums, and interpolation. Direct A) methods can be made as accurate as needed at the cost of increased computation, and serve as a benchmark in terms of accuracy. Typical examples of direct methods are described in P. M. Joseph, "An Improved Algorithm for Reprojecting Rays Through Pixel Images," IEEE Trans. Med. Imag. vol. 1, No. 3, pp. 192-198, 1982, D. C. Yu and S. C. Huang, "Study of Reprojection Methods in Terms of Their Resolution Loss and Sampling Errors," IEEE Trans. Nuc. Sci., vol. 40, No. 4, p 1174, August 1993, and G. L. Zeng, Y. L. Hsieh, and G. T. Gullberg, "A Rotating and Warping Projector/Backprojector for Fan-beam and Cone-beam Iterative Algorithm," IEEE Trans. Nuc. Sci., vol. 41, No. 6, pp. 2807-11, December 1994. However, in spite of various direct methods improvements, e.g., U.S. Pat. Nos. 5,559,335 to Zeng et al. and 5,625,190 to Crandall, direct methods require a high computational cost, proportional to N.sup.3 (denoted O(N.sup.3)), to generate N projections for an image with N.times.N pixels. For example, using one of these methods to compute the reprojection of a typical 4096.times.4096-pixel image requires 16.sup.3 =4096 times the computation needed for a 256.times.256 pixel image. Thus, various efforts have increased in academia and industry in an attempt to develop fast reprojection methods.
Another typical approach to accelerating the reprojection process in commercial products is to develop special purpose hardware that attempts to increase the rate at which the reprojection is performed. Such methods are described in E. B. Hinkle, J. L. C. Sanz, A. K. Jain, and D. Petkovic, "P/sup 3/E: New Life for Projection-based Image Processing," J. Parallel & Distrib. Comput., Vol. 4, No. 1, pp. 45-78, February 1987, T. M. Peters, "Algorithms for Fast Back- and Re-projection in Computed Tomography," IEEE Trans Nuc. Sci., vol. NS-28, No. 4, pp. 3641-3647, August 1981, J. L. C. Sanz and E. B. Hinkle, "Computing Projections of Digital Images in Image Processing Pipeline Architectures," IEEE Trans. Acoust., Speech & Signal Process., vol. ASSP-35, No. 2, pp. 198-207, February 1987, and U.S. Pat. Nos. 4,930,076 to Meckley, 5,008,822 to Brunnett el al., 5,136,660 to Flickner et al., 5,224,037 to Jones et al., and 5,559,335 to Zeng et al. Another approach described in C. R. Crawford, "Reprojection Using a Parallel Backprojector," Medical Physics, vol. 13, No. 4, pp. 480-483, 1986, and in U.S. Pat. No. 4,626,991 to Crawford et al., is to use a hardware backprojector, which is sometimes included in the reconstruction system, in a mode allowing the backprojector to perform a reprojection.
Such methods using special purpose hardware promise faster reconstructions than are obtainable using general purpose computers, but the special purpose hardware is expensive and offers only small factors of acceleration over less exotic implementations. The cost of implementing these methods, which typically split the work between multiple processing elements, is roughly proportional to the speedup they offer. Furthermore, with the astounding rate for increases in the performance of general purpose computers, the required special purpose architectures quickly become obsolete. For example, many modern systems no longer include a hardware backprojector, making the aforementioned method obsolete for these systems.
Another method, an analog optical processing method described in T. Lu, et al., "Projection Iterative Reconstruction Technique and Optoelectronic Implementation." Proc. IEEE ISCAS '92, San Diego, May 10-13, pp. 2469-2472, 1992, and U.S. Pat. No. 5,654,820 to Lu et al., uses analog optical signal processing to implement, among other functions, a tomographic reprojection. While this method may offer a large speedup, the method suffers from the same drawbacks listed above for other hardware acceleration schemes. Furthermore, being an analog processing technique, it includes well known limitations of optical analog signal processing such as limited dynamic range, nonlinearity, insufficient accuracy, and long-term drift. These limitations have restricted the use of analog optical signal processing, which has now been, for the most part, displaced by digital implementations in practically all fields of signal processing.
Other methods, such as Fast Fourier Transform (FFT) based methods proposed in C. R. Crawford, J. G. Colsher, N. J. Pelc, and A. H. R. Lonn, "High Speed Reprojection and its Application," Proc. SPIE--Int. Soc. Opt. Eng. Conf. Medical Imaging II, Newport Beach, Calif., vol 914, pt A, pp. 311-18, 1988, and U.S. Pat. No. 4,616,318 to Crawford, theoretically enjoy acceptable computational requirements of O(N.sup.2 log N). These methods, however, are based on a Fourier Slice-Projection Theorem, which is the basis for Fourier Reconstruction Algorithms (FRAs). As such, they suffer from the same deficiencies that have prevented FRAs from becoming commercially viable, although they have been known and rediscovered many times for more than 30 years now.
The Fourier Slice-Projection Theorem states that the Fourier Transform along the radial coordinate of a projection at angle .theta. is equivalent to a slice of the two dimensional Fourier transform of the image along the direction .theta.. A primary problem with this algorithm is the required step of interpolation between the rectangular grid in Fourier space on which the transform of the image is computed using the FFT, and the polar grid on which the Fourier transform of the projections must be evaluated. Poor theoretical properties manifest themselves in terms of severe, unacceptable artifacts in the reconstruction. The high amount of effort necessary to overcome these artifacts results in marginal or at best small, e.g., factor of 2-10, speedups over the direct techniques.
Another known method, a FFT-based reprojector-backprojector is a process for combined reprojection-backprojection, which is useful in some, but not all, iterative reconstruction algorithms. This method is described in A. H. Delaney and Y. Bresler, "A Fast and Accurate Iterative Reconstruction Algorithm for Parallel-beam Tomography," IEEE Trans. Image Process., pp. 740-753, May 1996. The method is based on a principle different from the FRA, but still uses the FFT. It is accurate, and its computational requirements are O(N.sup.2 log N). However, this algorithm is inherently restricted to a parallel-beam geometry and cannot be generalized to a fanbeam or other geometry prevalent in commercial diagnostic systems. Furthermore, the method cannot perform the reprojection operation in isolation, and therefore cannot perform many of the important functions of a reprojection algorithm.
Yet another method known as fast discrete Radon transform is a technique that performs a process similar, but not identical to reprojection. The method is described in M. L. Brady, "A Fast Discrete Approximation Algorithm for the Radon Transform," SIAM J. Comput., vol. 27, No. 1, pp. 107-19, February 1998. The method computes the discrete Radon transform. The method assumes a square pixel basis for the image, and only computes projections at a fixed set of non-uniformly spaced angles. A problem exists since the transform computed by this method is not a Radon transform as considered in general tomographic applications. Rather, it is defined in terms of partial sums of pixels whose center lies within a strip of predetermined width. As such, this method may have application in certain image processing algorithms, such as a Hough transform of binary images, but is inapplicable to tomographic reconstruction.
Still another method, reduced data processing as highlighted in U.S. Pat. No. 4,714,997 to Crawford et al., does not attempt to speed up the reprojection process itself. Instead, the approach reduces processing time by using an existing reprojection system with pre- and post-processing subsystems added to enhance the operation of the reprojector. In particular, this method operates on an image with reduced spatial resolution, i.e., fewer pixels, and possibly reduced size to reproject at fewer paths than the original data set. This method may provide adequate performance in certain artifact reduction applications utilizing reprojections, however, the speedup gains the method provides without compromising accuracy are small or even negligible for the high accuracy required of state-of-the-art systems.
Referring to FIG. 1, to demonstrate the need for speedup in known reprojection techniques, a square domain is shown with side lengths D. The square domain is defined such that a two-dimensional image .function..sub.c (x, y) vanishes outside the square domain, i.e., .vertline.x.vertline..ltoreq.D/2, .vertline.y.vertline..ltoreq.D/2. The image is approximated with a sufficient accuracy through interpolation of the image's N.times.N pixel discrete version .function.(m, n), i.e., ##EQU1##
where p is a pixel basis function and .DELTA.=D/N is the discretization size.
It is presumed that, in both x and y, the support of p is proportional to .DELTA. in size (denoted O(.DELTA.)), and p's essential bandwidth B is of O(.DELTA..sup.-1). Since p is spatially limited, p can only be essentially bandlimited, meaning that p has negligible energy outside the band. For example, the Fourier transform of an indicator function on a size .DELTA. square pixel has a main-lobe half-width of .DELTA..sup.-1 along each frequency axis, and approximately bandwidth B=.DELTA..sup.-1. Likewise, a tensor product of cubic B-splines used in numerical experiments discussed below has support 4.DELTA..times.4.DELTA. and bandwidth .DELTA..sup.-1. Alternative choices of essentially bandlimited pixel basis functions for tomography are discussed in R. M. Lewitt, "Multidimensional digital image representations using generalized Kaiser-Bessel window functions," J. Opt. Soc. Am. A, vol. 7, No. 10, pp. 1834-1846, October 1990. In practice, the standard pixel basis is often fi adequate, because the bandwidth of tomographic reconstruction is restricted by the filtering or convolution step of FBP.
A parallel-ray projection situation is depicted in FIG. 1. The Radon transform of the image EQU .function..sub.c (r,.theta.)=R.function..sub.c (r,.theta.)=.intg..function..sub.c (r cos .theta.-t sin .theta.,r sin .theta.+t cos .theta.)dt
produces parallel-ray projections of the image over a continuous range of angles 0.ltoreq..theta..ltoreq..pi.. This set of projections is often called a sinogram.
In many practical situations, the projections of .function..sub.c are measured only at a set of P angles {.theta..sub.i }.sup.P.sub.i=1, sometimes within a limited angular sector .vertline..theta..vertline..ltoreq..THETA., and O(N) uniformly spaced radial points {r.sub.k }. The required reprojection operation for a given image corresponds to the calculation specified by equation (2) for all {.theta..sub.i }.sup.P.sub.i=1 and {r.sub.k }.
Let EQU p(r,.theta.)=Rp (r,.theta.)
denote the projection of the pixel basis function. Clearly, it too has support of size O(.DELTA.). Using the representation of equation (1), it therefore follows that computation of .function.(r.sub.k, .theta..sub.i) at one view angle and one radial position requires O(N) operations. Thus, the operation count of the computation of projections at P view angles and N radial positions is of O(PN.sup.2), and for a typical case where P=O(N), the computational cost is of O(N.sup.3). This unfavorable cost scaling limits the use of the algorithm for many applications, and is unacceptable in iterative image reconstruction, e.g., reconstruction from limited data, which requires repeated reprojection of successive approximations of the image for angles {.theta..sub.i } for which the data is collected.
In view of the above, an object of the present invention is to provide an improved method that enjoys an accuracy of known direct methods, but requires less computational processing costs, i.e., a scaling of N.sup.2 log N or better. For typical image sizes this can result in speedup by a factor of 100 or more for a standard serial processor.
Another object of the present invention is to provide a cost-effective method that can use, but does not require special purpose hardware, and is easily implemented on standard serial and parallel architectures.
Yet another object of the present invention is to provide a method that enjoys speedup without the artifacts associated with methods such as FFT-based methods.
Still another object of the present invention is to provide a method that is flexible and can be adapted to the fanbeam or other geometries prevalent in diagnostic systems.